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We analyze the quantum melting of two-dimensional Wigner molecules (WM) in confined geome¬ 
tries with distinct symmetries and compare it with corresponding thermal melting. Our findings 
unfold complementary mechanisms that drive the quantum and thermal crossovers in a WM and 
show that the symmetry of the confinement plays no significant role in determining the quantum 
crossover scale nx ■ This is because the zero-point motion screens the boundary effects within short 
distances. The phase diagram as a function of thermal and quantum fluctuations determined from 
independent criteria is unique, and shows ‘melting’ from the WM to both the classical and quantum 
“liquids”. An intriguing signature of weakening liquidity with increasing temperature, T, is found 
in the extreme quantum regime. The crossover is associated with production of defects. However, 
these defects appear to play distinct roles in driving the quantum and thermal ‘melting’. Our study 
will help comprehending melting in a variety of experimental traps - from quantum dots to complex 
plasma. 
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I. INTRODUCTION 

A condensed phase of matter, and its transition to a differ¬ 
ent one is sharply defined only in the thermodynamic limitii^. 

Yet, the mechanism driving the transition and the “phases” 
can be discerned as a “crossover”^— even in a conhned ge¬ 
ometry with finite number of particles. Such has been the 
scenario for the “melting” of a Wigner molecule (WM)^ - 
a puddle consisting of Coulomb-interacting particles (elec¬ 
trons). Under the influence of thermal or quantum fluctu¬ 
ations, a WM crosses over to a liquid-like phase. WM is 
a nano-scale version of two-dimensional (2D) Wigner crys- 
tal^ii^, an insulating phase of a 2D electron solid whose non¬ 
conducting nature arises from the propensity of the electrons 
to localize as far apart from each other as possible (consistent 
with their density) minimizing the free energy. 

The study of the Wigner molecules found a large scientific 
quest primarily with the advances in nano-technology. They 
are easily tunable using electrostatic and magnetic meth- 
odsii employing non-invasive techniques. Their importance 
in fundamental Physics is crucial, providing a hotbed for 
studying the complex interplay of Coulomb-repulsion, quan¬ 
tum interference effects in the confinement, level quantiza¬ 
tion due to their smallness, and finally, the disorder in the 
form of irregularities in the confinement geometry. 

The recent proposal of a insulator-metal transition in 2D 
electron gas (with inherent disorder)^ has renewed the in¬ 
terests for the study of “disordered Wigner melting”— lii, 
because such transition is often attributedi^ to a melting of 
a Wigner crystal or a Wigner glass. Disorder or impurities 


break all the spatial symmetries of the pure system reducing 
the ability of the particles to delocalize. There have been pro¬ 
posals for new and intervening phasesi^^— between the crys¬ 
tal and liquid. Finally, the presence of disorder can change 
the nature and criticality of a transitio n^^i^° , and such mod¬ 
ifications of Wigner melting have not yet been looked into. 
Overall, the complex interplay associated with such melting 
remains as one of the outstanding problems of condensed 
matter physics. 

Seeking the physics of Wigner melting in confined geome¬ 
tries has significant experimental relevance^! in varied areas, 
such as, radio-frequency ion traps^^, electrons on the surface 
of liquid He^^, electrons in quantum dots in semiconductor 
heterostructures^l^ and in dusty plasmas^i. Unlike in the 
bulk system, the long-range Coulomb (~ r~^) repulsion is 
poorly screened in a finite cluster. Experimental clusters, 
such as large lateral quantum dots2^ often have inherent ir¬ 
regularities of the confinement, which are expected to act as 
disorder. This conjecture found support from the tunneling 
conductance measurements of chaotic dots^, and associated 
Coulomb blockade experiment a^^i^^ . Although a quantum 
melting of Wigner-type in a circular trap has been studied ex- 
tensivelj*^!^'^, a similar study in the irregular confinement 
has not yet received as much attention, see however Ref—. 
Recently, two of us reported the classical “melting” of irreg¬ 
ular Wigner molecule (IWM) under the influence of thermal 
fluctuations^, devised criteria to quantify such crossover, 
and identihed its mechanism in terms of defects. There are 
fundamental questions to address upon inclusion of quantum 
fluctuations in the form of zero-point motion. How different 
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is the nature of such melting in a IWM from that in confine¬ 
ments with certain symmetries? What roles do the defects 
play for such a quantum crossover? Finally, what would be 
the phase diagram in the plane of thermal and quantum fluc¬ 
tuations? Addressing such questions in the present paper, we 
list below our main findings: 

(1) The symmetry of the confinement does not control the 
“solid” to “liquid” crossover scale in WM, neither does it 
influence the underlying “melting” mechanism. 

(2) The crossover turns out to be unique irrespective of the 
used distinct criteria for “melting”. 

(3) The phase diagram consists of three phases: a classical 
liquid, Wigner molecule and a quantum liquid along with in¬ 
triguing re-entrant solidification with increasing temperature 
in the extreme quantum regime. 

(4) The crossover in the traps with spatial symmetries indi¬ 
cates that in the thermodynamic limit the quantum melting 
is likely to be of first order with the characteristic disconti¬ 
nuity in the observables. The weaker melting in an irregu¬ 
lar confinement is more subtle for predicting thermodynamic 
limit. The thermal melting, on the other hand, follows the 
disclination mediated KTHNY mechanism^^r^. 

The rest of the paper is organized as follows. In Section II 
we introduce the model, the parameter space and the sim¬ 
ulation method used to study both quantum and classical 
systems. In Section III, we present our main results. The 
corresponding subsections include the comparison of several 
melting criteria, the reconstructed diagram and the discus¬ 
sion of the melting mechanism. 


II. MODEL, PARAMETERS AND METHODS 

Our model of WM consists of N distinguishable particles 
(Boltzmannons) in 2D plane and interacting via long-range 
Coulomb repulsion. These particles are trapped by an exter¬ 
nal potential I4onf(r)- For much of our calculations we com¬ 
pare results from three confinements: (a) Circular, V^onf(r), 
(b) Elliptic, I4E;jf(r) and (c) Irregular, V;^o„f(r), given by, 

Konf(j') = ar^ 

Konf(r) = a{x'^/b + by‘^-2Xx^y^+'y{x-y)xyr} (I) 

where r = \Jx^ y"^ and a = imwg. The confinement 
strengths a, a and c, d controls the average particle density in 
these traps. Because we expect the (’’) mimic the uni¬ 
versal features of disordered systems, it is designed to break 
all the spatial symmetries. For example, b breaks the x-y 
symmetry while A introduces chaoticity and 7 breaks reflec¬ 
tion symmetry. Appropriate values for parameters represent¬ 
ing universal disorder physics are found in literatur o^°i^^~ — . 
While disorder is typically introduced in bulk systems by 


random impurities, they originate in a confined system with 
small number of particles from the irregularities at the ‘soft’ 
boundary which are relevant for experiments in which the 
external confinement is set up by electrostatic and magnetic 
meansii. The ensemble of particles in this potential will be 
referenced below as IWM. 

In case of a harmonic trap, the total Hamiltonian is given 

by 
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where m and e are the effective mass and the dielectric 
constant of the medium, respectively. Next, we introduce 
the length and energy scales rg and Eg, specified by 
e^/erg = mL>j^r'^l2 and Eg = e^/erg. After the scaling 
transformation {r —^ r/rg and E —>■ E/Eg} Eq. ([2|) takes 
the dimensionless forni^ 
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Here n = \f2l\lr'^ is the quantum parameter related with 
the strength of the trap defined via the oscillator length l^ = 
h/mojQ. The analogy to the bulk in deciding the quantum 
parameter is given by where rs specifies the 

average electron density, see Eq. ©■ The thermodynamic 
properties of ([3]) have been studied in the canonical ensemble 
(hxed particle number N and temperature T = ksT/Ef). 

The potentials and are quadratic in r, 

whereas, the V/onf(’') bas quartic dependence. To match the 
dimensions of three potentials, we rescale the strength of 
the irregular potential as a —>■ a' (mw^/2rg). This rescaling 
ensures that {a', b, A, 7 } in Eq. (1) can be treated as dimen¬ 
sionless parameters, whereas the definition of the quantum 
parameter n, being the pre-factor in the kinetic energy term, 
remains the same for all the three confinements. 

The thermal melting of an irregular Wigner molecule has 
already been studied in detail^. In the current study we 
focus on the effect of quantum fluctuations (by varying n). 
Temperature T is kept fixed at a low value, corresponding to 
the solid phase in the classical regime. By increasing n we 
induce quantum melting in all the three confinements. 

The quantum A^-particle system ([3]) can be analyzed by 
first principle simulations using the path integral Monte 
Carlo (PIMC) techniqu o^^d° . The particle statistics is lim¬ 
ited to “boltzmannons” (no exchange) to distinguish the ef¬ 
fects related with the disorder and quantum statistics. In 
the path integral representation the A-particle density ma¬ 
trix is mapped on that of a classical system of interacting 
“polymers”. Particles are represented by the trajectories 
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which evolve in the imaginary time 0 < t < /3, with the up¬ 
per bound specified by the inverse temperature /3 = l/fcsT. 
Number of time-slices M changes with /3 and the quantum 
parameter n by the relation M = ipn, where I is an integer 
in the range 1 to 10 to reach convergence in the energy better 
than 5%. In most of our simulations we use M = 100 inde¬ 
pendent of temperature T. The ensemble averages have been 
estimated over ~ 10 ® independent realizations of particle tra¬ 
jectories efficiently sampled via the bisection algorithm's. 

Complementary, the classical counterpart has been stud¬ 
ied by standard simulated annealing Monte Carlo (MC) 
schemeSi based on the Metropolis algorithmSS. In the simu¬ 
lated annealing schedule we start from higher T and anneal 
down to the desired T to track the appropriate low energy 
states. 

It is important to note, that the identification of the melt¬ 
ing transition in the case of irregular trapsSi filled with 
quantum particles is more complicated compared with the 
classical clusters formed in circular symmetric potentials. 
The latter have been successfully studied via the modihed 
Lindemann-based^ parameters, the orientational and radial 
correlation functions which take into account coexistence 
of two different types of symmetries, i.e. ordering into a 
triangular-lattice structure (with hexagonal symmetry) in 
the bulk of the cluster and a shell-structure in the boundary 
region^^— , as wells as possible radial inhomogeneity of the 
onset of melting. 

In contrast, the onset of melting in quantum clusters will 
be signihcantly masked by the zero-point fluctuations. In 
particular, the Lindemann-based parameters and the corre¬ 
lation functions in quantum solids at T = 0 get a significant 
background shift related with the quantum mechanical un¬ 
certainty in the particle position. As an example, the critical 
value of the Lindemann parameter at the solid-liquid quan¬ 
tum transition^ (T = 0 ) is increased to 0.2 from its origi¬ 
nal value O.I valid for classical solids. Similar observations 
hold for finite quantum clusters^i^i^. The combined effect 
of quantum fluctuations and the order-disorder transition is 
even more difficult to distinguish in small systems {N < 100) 
where a sharp transition is replaced by a crossover behavior. 

Partially to overcome this problem in the present studies, 
we measure the order parameters characterizing the “solid 
phase” in two different ways. First, by taking the ensemble 
average over particle configurations on every imaginary time 
slice. This is called head-by-bead (BB) distribution and corre¬ 
sponds to the correct quantum mechanical average. Second, 
we exclude the zero-point fluctuations related with the delo¬ 
calization of quantum particles and address the bead-centroid 
(BC) positions during the imaginary time evolution 

I 

(r*) = ^ / r,(r)dr, i = l,N. (4) 

P Jo 

The physical picture is captured the BB-distribution as 


it takes into account the wave function overlap. The BC- 
distribution, instead, is more suited to describe different 
structural symmetries, cf. relative particle positions on the 
shells, and corresponds to the semi-classical picture. The 
BC-picture also depicts the evolution of particles in the con¬ 
figuration space, while the BB the combined evolution in¬ 
cluding the imaginary time. Because BC-calculations pro¬ 
vide a semi-classical interpretation for the quantum effects, 
it is useful for the comparison between the classical and its 
corresponding quantum counterpart. Both pictures are nec¬ 
essary to fully understand different melting criteria intro¬ 
duced below. 


III. RESULTS 

In our studies we mainly concentrate on the system with 
N = 57 particles. In the case of circular potential it is called 
Magic Cluster as the hexagonal symmetry dominates over 
the full structure^. In irregular and elliptical case, no such 
concept as a magic cluster exists. 

Next, we should address the question: which parameter 
needs to be fixed to compare melting in three different traps 
on equal grounds? The parameter n is inversely proportional 
to ro, the spacing between two nearest neighbors. We tune 
the trap parameters (c and d in the elliptical trap, and a' in 
the irregular one) to make rp equal in the three cases. This 
results in the same average density. To do this, we com¬ 
pare the coordinate rp for the ground state of two particles 
obtained by minimizing the potential for the different traps. 

The trap parameters found by this procedure are: a' = 
0.0954 in case of irregular, and c = 1.25 and d = 0.85 in case 
of elliptical potential. The other parameters for the irregular 
trap are similar to those used in Refi^, cf. & = 7 r/ 4 , A = 
0.635 and 7 = 0.2. How close the spacing rp approximates 
the average inter-particle distance can be checked from the 
pair distribution functions presented below. 


A. Study of the crossover from density profiles 

Clear evidence on the melting transition in finite systems 
(“solid”-“liquid” crossover) can be directly gained from two- 
dimensional density distributions. Fig. [1] presents a sequence 
of density plots for the IWM in a given confinement. 

The left column presents an example of melting by ther¬ 
mal fluctuations (the quantum parameter has the smallest 
value n = 0.01). Here, fluctuations are nearly frozen at low 
T and n and individual particles are resolved as bright spots 
in the density profile (though certain particles show quantum 
signature through a spread in their wave functions). With 
the increase in T, particle delocalization occurs in two ways: 
(a) Thermal diffusion of particles around their equilibrium 
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n=0.01 n=0.06 n=0.10 n=0.15 



FIG. 1: 2D density p(x,y) plots on the T-n plane (temperature- 
quantum parameter). Quantum fluctuations n are increased from 
left (n = 0.01) to right (n = 0.15) column. Following along a 
column (fixed n) or a row (fixed T) one can trace the melting in¬ 
duced either by thermal or quantum fluctuations. Note, that due 
to the combined effect of the external field and pair interactions 
the density fluctuations are strongly anisotropic. A typical cluster 
size L/ro can be read out from Fig. [O] Our length unit scales as 
ro oc 1/n^. Thus, while the rescaled system size stays practically 
unchanged with n, its absolute spatial dimensions evolve quite 
rapidly with n. 


position, and (b) The thermal fluctuations create occasional 
long and string-like paths of delocalization, whose role on 
the crossover is already discussed in the literature^. Mov¬ 
ing along a row, (say, the one corresponding to T = 0.005), 
the strength of the quantum fluctuation increases and the 
corresponding crossover occurs primarily by the spreading of 
wave functions of the individual particles. We found no ad¬ 
ditional ‘string-like paths’ upon increasing n. However, such 
paths, once developed for T > 0.01, survive quantum fluc¬ 
tuations up to a large n and the melting becomes complete 
when such paths cannot be resolved from the spread of par¬ 
ticles’ wave function. This is clearly seen in the third row 
from the top, corresponding to T = 0.020. 

Now let us draw attention to the fV-particle wave function 


(or density) at n = 0.15 and T = 0.005. Here the individ¬ 
ual particle positions are hardly distinguished resulting in a 
quantum liquid state. Interestingly, a rise in the temperature 
to r = 0.010 leads to slight but noticeable particle localiza¬ 
tion in the center and the first row. This illustrates a gradual 
transition from quantum to quasi-classical picture. Maxi¬ 
mum of both thermal and quantum fluctuations is reached 
in the bottom right panel (in our parameter space) and rep¬ 
resents a ‘melted’ IWM. 

While the qualitative picture of the crossover is appar¬ 
ent from Fig. [TJ few points deserve special mention: (a) In 
the quantum cluster at n = 0.15 the particles are reason¬ 
ably delocalized even at the lowest T. However, heating to 
T = 0.020 seems to weaken the “liquidity” and distinction 
of individual particles becomes better than for T = 0.005. 
This is in contrast with the scenario common for low n 
(n = 0.01), where the temperature increase ‘dissolves’ the 
particles monotonically. This observation indicates emer¬ 
gence of the incipient re-entrant behavior showing weakening 
“liquidity” with increasing T. A final classical liquid is ob¬ 
viously expected for T —>■ oo. This feature is qualitatively 
consistent with the similar trend observed in bulk systems, 
where re-entrant crystallization from a quantum liquid phase 
by increasing the temperature was reported^; (b) As can be 
clearly seen from Fig. 1 all particles do not melt simulta¬ 
neously upon increasing the fluctuations. The particles on 
the ‘string-like path’ diffuse in the background more easily, 
because, the spreading of wave function is easier as these 
particles undergo larger thermal displacements; 

(c) In our rescaled units the average inter-particle distance 
stays practically unchanged with variation of n, i.e. (r) /ro ~ 
1. This also holds for the cluster size, see Fig. [T] However, 
in the absolute units (r) shows a strong n-dependence: (r) ~ 
ro = as/n^. Similar dependence is captured by the quantum 
Bruckner parameter 

r* = = (r) /[osv^ = l/[y/TTn\ (5) 

typically used to characterize a bulk density of homogeneous 
systems. Finally, Eq. (O can be used for comparison of the 
onset of melting in bulk and trapped systems. 


B. Crossover in bond orientational order 

The crossover in finite systems bringing out the correlation 
driven Wigner physics, naturally warrants proper qualifica¬ 
tion. Here we study the evolution of the bond orientational 
order (BOO) ^°i^^ with n. The hexagonal symmetry and de¬ 
viation from it can be quantified by the order parameter 

= (V’eW) = X! exp(6z6»„„(r))\ , (6) 

\ nn—1 / 








































5 



FIG. 2: (a-c) Distribution of the bond-angular order parameter 

P{ip 6 ) shown for different strength of quantum fluctuations n. vE'e 
is evaluated from the bead-centroid (BC) particle coordinates (ri), 
Eq.®. Three different traps (a)-(c) are analyzed at hxed temper¬ 
ature T = 0.005. Common behavior is observed: the maximum, 
Pn ~ max[P('!/)6)] (n = 1,2), gradually shifts from the position 
Te ~ 0.9 specifying a “solid” phase to Te ~ 0.3 typical for a 
quantum or classical “liquid”. The weakness of bond-orientation 
order in the irregular trap reflects in the weak peak at Pi even 
at the lowest n, and P{tpe) in an irregular trap features a second 
peak drifting towards P2 with increasing n. d) The n-dependence 
of Pd{n) = Pi — P2 during the “solid”-“liquid” transition. Pd{n) 
demonstrates a similar trend for all three traps and changes its 
sign around tie ~ 0.10(1). The inset shows the evolution of Pd in 
a similar classical systems. Here it is concave in nature while in 
the quantum case it is of convex nature. 

where the sum is taken over six nearest neighbors [nn) of 
a particle located at the position r, and is the relative 
angle between the vector — r and an arbitrary fixed axis. 
The ensemble average is taken over all particle positions. 

Evaluation of 'I'e has been performed using the BB- and 
BC-methods (see above). In the first case, the BOO includes 
quantum mechanical fluctuations, and as a result ffg is sig¬ 
nificantly reduced in a quantum solid even with a perfect 
hexagonal symmetry. In the bead-centroid (BC) method the 
averaged particle positions are used, allowing a semi- 


classical interpretation and demonstrating the behavior typ¬ 
ical for classical solids. 

Fig. [ 2 ] compares the distribution of the order parameter 
'I'e = 'I'g^ for three different confinement potentials, several 
values of n and lowest temperature T = 0.005. In the en¬ 
semble average ([H]) we excluded the particles at the cluster 
boundary. 

A comparative study of P(ll'6) in the three confinements 
brings out their common qualitative evolution: At small 
value of quantum parameter (n = 0.01) ^(^'e) is peaked 
around ife = 0.9 indicating that in the cluster bulk the ma¬ 
jority of particles have six nn with the hexagonal symmetry. 
As n is increased the peak weakens and shifts progressively to 
lower values. Finally, at n = 0.15 we observe a rather broad 
distribution centered around dig = 0.3. Further increase of n 
does not significantly alter this shape. Inspite of these broad 
similarities, the evolution of in the irregular trap dif¬ 

fers from those in a circular or elliptic confinement: Each 
trace of features double peaks up to n ~ 0.10 and the 

one at ikg = 0-9 is weaker than those in the traps with spa¬ 
tial symmetries. The second peak, far weaker than the first 
one, drifts towards ike = 0.3 contributing to the sole peak at 
large n. 

The convergence of the distribution P{'i>e) to one charac¬ 
teristic shape for all confining geometries at large n can be 
explained by the fact that any peculiarities of an external 
potential are effectively screened due to quantum delocaliza¬ 
tion of particles. This induces the effect of smoothness in 
the external field as has been demonstrated by the varia¬ 
tional perturbation theory due to Feynman and Kleinert^. 
The part of a full partition function related with an exter¬ 
nal field can be substituted by an effective classical potential 
which accounts for possible effects of quantum fluctuations. 
The induced smoothness, also in the pairwise interactions, is 
well known in the literature^. 

In the limit of small n we also observe quite universal form 
of P(ll'6) with the peak at d'g = 0.9. This again comes from 
the screening effect but classical in its origin. It is specific 
to Coulomb systems, that particles mainly experience inter¬ 
action with their nearest neighbors, whereas the interactions 
at larger distances coming from the outer parts of the system 
are canceled. The cancellation will be complete in the mean 
field approximation. 

Another common feature found for all three traps is the bi- 
modal nature of the distribution. The height of ^(d'g) at two 
reference values, 'kg = 0.9 and ikg = 0.3, which distinguish 
“solid” and “liquid” phases, can be used to characterize the 
degree of “solidity” and “liquidity” in the crossover regime. 
We have evaluated the corresponding heights Pi = P(0.9) 
and P 2 = P(0.3) and their difference Pd for 0.01 < n < 0.15. 
We plot Pd in Fig.[2]for all three types of the confinement and 
observe a qualitatively similar behavior. The n-dependence 
of Pd demonstrates a monotonic decay from positive to neg- 
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ative values. This allows to identify the crossover from a 
“solid” to a “liquid” phase. The point nx ~ 0.10(1) [which 
corresponds to the bulk density ^ 56(10) using Eq. JS])] 
where Pd changes the sign is the mid-point of the crossover. 
Similar to our previous discussion of the BOO, the nature 
of the confinement does not play a significant role in the 
transition induced by quantum fluctuations. 

The inset in the Fig. [2] presents similar analyses for the 
classical system in the same confinement (of all three types) 
where the transition is triggered by thermal fluctuations. 
The comparison shows two distinct differences. First, al¬ 
though the evolution of Pd with n (at fixed T) is of convex 
nature, Pd versus T (at n = 0) has a concave curvature. 
Second, the BOO in V[.Qnf(r) is destroyed earlier compared 
to V[,Q[jf(r) and V[,Q[jf(r) . Both these features can be un¬ 
derstood by connecting the mechanism of melting with the 
defect production, which will be discussed in the details in 
Section G. 

For completeness we also calculated ^(^'e) from dH]) us¬ 
ing the BB-distribution, dig = dig®. In this approach we 
capture additional (background) fluctuations from the indi¬ 
vidual particle trajectories in the imaginary time. Although 
P(di'g®) shows similar trend as T’(d'f ^) its larger background 
fluctuations makes the identification of the crossover more 
difficult. 

It is interesting to check whether the re-entrant behavior, 
as observed from the density plots in Fig. [I] can be more 
quantitatively analyzed via the BOO. In Fig. [3] we present 
the T-dependence of P('I'g) at the largest value of quantum 
parameter n. In this case, surprisingly, the “liquidity” of 
the system defined in terms of the height of the distribu¬ 
tion P 2 = T’(d' 6 )|® 6 =o .3 dominates at low temperatures. As 
the temperature is increased from T = 0.005 to T = 0.04 
the height P 2 decreases, whereas the right wing of the dis¬ 
tribution around d'g = 0.9, characterizing the degree of “so¬ 
lidity” (or hexagonal symmetry), steadily increases. This 
leads to the conclusion that the hexagonal order is higher at 
T = 0.04 then T = 0.005. This finding can be explained by 
the anisotropy of quantum fluctuations at low temperatures. 
Quantum particles are mainly delocalized along narrow po¬ 
tential valleys, following the shape of the trap, as shows the 
comparison of the density plots in Fig. [1] for n = 0.10 and 
n = 0.15 at T = 0.005. This naturally leads to a significant 
distortion of the BOO as we approach the cluster bound¬ 
ary, where the symmetry of the irregular potential dominates 
over the hexagonal order. In contrast, at high temperatures 
the particles get more localized and the BOO is dominated 
by the bead-centroid particle coordinates (4) which resemble 
the hexagonal order with some distortions. As a net result 
of this transition to the quasi-classical regime the hexagonal 
symmetry is partially restored. 

The lower inset in Fig. |3] shows the difference of the peak 
heights Pd of two “phases”. Starting from negative values. 



FIG. 3: Main panel: T-dependence of P(v['6) deep in the quan¬ 
tum regime (n = 0.15). The re-entrant behavior is clearly ob¬ 
served for 0.01 < T < 0.02. Lower inset: T-dependence of the 
height difference Pd- The non-monotonic behavior is observed in 
the same temperature range and is related with a partial restora¬ 
tion of the hexagonal order. Upper inset: Typical behavior of 
P{tp6) in the quasi-classical case (n = 0.01). A similar depen¬ 
dence is observed in the classical limit (n = 0). 


due to a loss of the hexagonal order in the boundary region at 
low T, Pd shows a non-monotonic T-dependence. It increases 
for the temperature range 0.01 < T < 0.02 being a clear 
indication of the weakening of “liquidity” and resembling the 
re-entrant behavior observed in the macroscopic systems^. 
However, we should stress that in our case this feature is 
observed due to a strong perturbation of the BOO in the 
boundary region at low temperatures. As discussed above 
the BOO is restored at high temperatures. The prominence 
of this effect is naturally expected in finite size systems where 
the boundary effects play an important role. 

Finally, we analyze the evolution of the BOO in the quasi- 
classical regime where the thermal fluctuations make impor¬ 
tant contributions. The upper inset in Fig. |3] shows the T- 
dependence of Piipe) at n = 0.01. We find the standard 
behavior as observed in Fig. [2] during the crossover Pd de¬ 
creases monotonically and there is no sign of the re-entrant 













7 


behavior in ^(^g)- 


C. Study of the crossover via the Lindemann ratio 


An additional quantity which tracks the solid-liquid 
crossover is the Lindemann parameter— 
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In the ensemble average (...) we use the bead-by-bead dis¬ 
tribution of particle trajectories from the PIMC simulations. 
Here M is the total number of time slices and Vu is the posi¬ 
tion of the i-th particle at the t-th time slice relative to the 
trap center. 

In Fig. m^a) we show the evolution of with n at low¬ 
est temperature (T = 0.005) for all three traps. First, 
evolves linearly with n, and a distinct change in the slope 
takes place only around nx ~ 0.11(1). This holds for all 
three cases validating that the quantum crossover param¬ 
eter nx is universal for different trap geometries, though 
the change of slope in the irregular confinement is rounded 
off beyond nx- We note that demonstrates a strong 
dependence on the quantum fluctuations n even in a solid 
“phase”. This is in contrast to the findings for classical sys¬ 
tems, where only a weak dependence on thermal fluctuations 
was observed below a critical temperature (T < Tx) ^^'^^ - 
Our observation is not surprising, as in the quantum case 
the Lindemann parameter captures the background contri¬ 
bution coming from the evolution in the imaginary time and 
the quantum mechanical uncertainty. Both are growing with 
n. In fact, (not shown here) that picks up only a semi- 
classical contribution, shows more resemblance with the clas¬ 
sical behavior as expected. 



FIG. 4: (a) Evolution of with n for three different traps. 
In all cases we observe a similar trend with a characteristic kink 
which distinguishes two slopes in the evolution of with n, 
and defines nx{T). Such differences in the slope become weaker 
with T. (b) The locns of nx{T) traces the curve ABCDEF and 
translates into the phase diagram in the T — n plane shown in (c). 
Similarly the evolution of Pd with T at (d) gives more confidence 
to draw the phase diagram, (c) Phase diagram: two curves rep¬ 
resenting nx{T) (estimated from two independent criteria) are 
compared and distinguish different phases. 


D. Phase Diagram 

Our detailed analyses of versus n for different temper¬ 
atures, cf. Fig.IlKb), allows us to generate a phase diagram, 
cf. Fig. 4(c). The locus of nx{T) traces the curve ABCDEF 
in Fig. |4l(b) and is translated into the phase boundary in the 
T—n plane. 

To strengthen our confidence on the general validity of 
the phase diagrams reconstructed based on the Lindemann 
parameter—, we perform, in addition, an independent recon¬ 
struction using our Pd-criterion, cf. Fig. [21 The critical 
points nx(T) are defined by the sign change in P^. As an ex¬ 
ample, Fig. |4l(d) shows Pd versus n for several temperatures 
T in the case of the V),Qjjf(r). 

The phase boundaries reconstructed from both criteria 
demonstrate good agreement, cf. Fig. |4l(c). We can con¬ 


clude, that the use of both methods for the melting-analyses 
is justified and mutually complementary. 

Next, we discuss some important features of the phase 
diagram. The left-most side of the phase boundary (for 
n < 0.05) corresponds to a ‘quasi-classical’ regime where 
the degree of “liquidity” is mainly due to thermal fluctua¬ 
tions. In contrast, quantum fluctuations play a dominant 
role for the right-most phase boundary and induce the zero- 
temperature melting around nx ^ 0.11(1). 

The slope of the phase boundary for small n can be well 
described by a single coupling parameter F well known for 
classical systems. The tangent is given by F“^ = ksTd/iV), 
i.e. the classical melting temperature measured in units of 
the characteristic potential energy {V). Quantum fluctua¬ 
tions systematically reduce Pci, as the system deviates pro- 













gressively from the classical limit, making the boundary line 
sub-linear in n. Deep in the quantum regime, when the zero- 
point fluctuations dominate, the curve bends down fast. The 
maximum in the melting temperature is reached for some in¬ 
termediate n-values. Remember, that n was defined as the 
reduced density parameter scaling as n oc 1/rg. Hence, the 
absolute value of the critical temperature shows a non¬ 
monotonic behavior. First, Tc is increased with the density. 
Then, above some critical density (n* ~ 0.09), it fast reduces 
to zero. Clearly, the observed peak in the reduced tempera¬ 
ture T = ksT / {V) around n* occurs at T < Tci{n*). 


E. Correlation functions 

We further quantify the solid-liquid crossover and the dif¬ 
ferent phases from the study of the spatial pair correlation 
function^ 
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and the bond-orientational correlation function 


ffeW = (V^edr'DV’edr' - r|)) • 


( 8 ) 

(9) 


The spatial decay of ge{r) characterizes the correlation 
length of local orientational order tjje{r), while g{r) describes 
the probability to find any pair of particles at distance r. 

The amorphous solids (as IWM) is characterized by its 
quasi-long range orientational order and a depleted posi¬ 
tional order— in the limit oi n,T —?► 0. It is interesting 
to study how these correlation functions evolve with n, i.e. 
in the presence of quantum fluctuations. Such results are 
presented in Fig. [51 The function g{r) was calculated by the 
BB-method, while ge{r) using the BC-method [in Eq. m we 
substitute '06 (r) = 0f®(r)]. 

Our results for ge{r) (left column in Fig. [5]) show well de¬ 
fined Bragg peaks up to a distances comparable with the 
linear dimension of the system, indicating the long-range na¬ 
ture of the BOO in all three traps. By increasing n, the 
Bragg peaks broaden progressively up to nx ~ 0.11. Then 
they disappear abruptly leaving bump-like features of a liq¬ 
uid. This is consistent with the crossover scale inferred from 
other criteria introduced in Sec. Hira and ImCl The long 
range nature of these correlations and their sudden demise 
are independent of the nature and symmetry of the confine¬ 
ment. 

For the circular and elliptic confinements g{r) provides 
sufficient information on the positional order. The Bragg 
peaks of such a mesoscopic “solid” phase can be well resolved. 
The disappearance of the peaks with quantum fluctuations 
is smooth and is in contrast to the behavior of g&{r). 



FIG. 5: Left column: Evolution of geir) versus n estimated with 
the BC-method for three different confinements. Right column: 
Evolution of the pair correlation function g{r). The legend indi¬ 
cates the n-values used in all plots. The center of the crossover is 
identified with nx = 0.10 when the Bragg peaks in gc{r) disap¬ 
pear. Similar observation holds for all three traps. 


Slightly different behavior is observed for the irregular 
trap. Here, the peaks in g{r) are broader and overlap even 
for small n-values indicating that there is no positional order. 

Finally, we note that the location of the first peak in gc,{r) 
and g{r) is close to one (in units of rg) irrespective of the 
confinement type. This validates our choice of parameters for 
the trapping potentials to ensure the same average density. 
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F. Signature of the classical/quantum crossover in the 
fraction of defects 

Thermal melting in 2D follows the KTHNY— theory 
and demonstrates some universal features mediated by the 
unbinding of topological defects. It is the dislocations and 
disclinations^ that signal the destruction of the positional 
and orientational orders, respectively. The irregularities in 
the confined systems, if any, are known to mask the prolifera¬ 
tion of the free dislocations^, but the disclinations still play 
a major role in the thermal melting. The natural question 
thus arises: What role do the disclinations play, if at all, for 
the quantum melting in confined geometries? 

We address this question by analyzing the fraction of par¬ 
ticles Nc^/N with the coordinate numbers Cn (n = 4,... 8) 
for a wide range of the quantum parameter n and compare 
the results with thermally induced melting. 

The particle numbers Nc^ have been unbiased estimated 
by the Voronoi construction^. The results for the three traps 
are presented in Fig. [6] quantum melting (upper row) and 
thermal melting (lower row). In both cases, Nc^ makes the 
largest contribution in the “solid” phase for all the traps. 
This fraction, however, is steadily decreased by both thermal 
and quantum fluctuations in the expense of generating new 
disclinations, i.e. increase of the populations Nc^ with Cn ^ 
Ce. 

The fractional change in the number of these defects for 
circular or elliptic conhnements are more abrupt as a func¬ 
tion of n (upper panel) compared to their gradual evolution 
with T (lower panel). The convex (or concave) nature of the 
evolution of Pd-, at least for the circular and elliptic traps 
in Fig. [2] in quantum (or classical) case can be understood 
by the nature of evolution of the defect-fraction. A simi¬ 
lar argument for the irregular trap is masked by a rather 
weak Pd, even deep in the “solid”, and makes such a cor¬ 
relation between the shape of Pd and Nce less substantial. 
The lack of symmetry in the V))onf(’') leads to a relatively 
fast destabilization of six-coordinated neighbor which makes 
Pd to approach zero at a slightly lower T and also a lower n 
than in other traps. Though this difference in the tuning pa¬ 
rameter for classical and quantum fluctuations stays within 
the crossover interval ATjc or Anx- The convex nature of 
Pd vs n as opposed to the concave nature of Pd vs T basically 
follows from the similar nature of the loss of six coordinated 
neighborhood. 

It will be erroneous, however, to conclude from Fig.IHlthat 
the crossover mechanism has a similarity for the thermal and 
quantum melting. For description of its critical properties a 
d-dimensional quantum system can be mapped onto the cor¬ 
responding (d -I- l)-dimensional effective classical system. A 
prominent example is the quantum Ising model. The size 
of the extra dimension is determined by the inverse temper¬ 
ature (3 = l/ksT and diverges as T —^ 0. Hence, at low 



FIG. 6: Fraction of the defects versus the tunning parameter: 
n in quantum (upper row) and T in classical (lower row) cases. 
Qualitatively similar trends are observed: the particle fraction, 
Nc„/N, with the coordinate numbers Cn (n = 5,7) steadily in¬ 
creases with both thermal and quantum fluctuations. The decay 
of the hexagonal order, cf. Nc^/N, is of convex nature in the 
quantum case though the curvature is rather weak for the irregu¬ 
lar confinement as well as its magnitude at the lowest n. In con¬ 
trast, its evolution is of concave nature in the classical case. The 
mechanism behind melting in the classical system is the prolifer¬ 
ation of disclinations. It is not conclusive for the quantum case 
due to the evolution in the imaginary time and effective presence 
of third dimension. 


temperatures quantum effects are always important. As a 
result, in a 2D macroscopic system the order of the phase 
transition is changed: second order for the thermal melting 
(KTHNY-type) is substituted by first order for the quantum 
melting (variation of the density parameter at fixed low 
r)iii^. Detailed calculations done for the quantum 2D XY 
model provide the following scenario^. The correlation func¬ 
tion shows the crossover between the 2D and 3D behaviors at 
Rc ~ 1/T. For the distances R> Rc a, system exhibits only 
the 2D character. Up to the distances R < Rc there exists 
a quasi-long-range order and the correlation function shows 
the 3D behavior. Once the system size satisfies L ^ Rc the 
3d character dominates. By lowering the temperature, Rc 
diverges and the 3D behavior can, in principle, be observed 
in an arbitrary system of finite size. 


Classical Quantum 
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FIG. 7: (color online) Particle trajectories are shown in dif¬ 
ferent confinements to contrast the mechanism of quantum and 
classical melting. The top two rows show the trajectory of the 
BC-coordinates of the particles in irregular and elliptic confine¬ 
ments for three different n values, whereas, the bottom two rows 
show the same during classical melting at three T. In the classical 
IWM, melting starts by coherent movement of certain particles 
along a tortuous path while other particles rattle around their 
equilibrium locations. In the quantum counterpart, we observe 
harmonic motion of all particles around their equilibrium posi¬ 
tion brings in the quantum melting. Colored lines show how the 
individual particles diffuse into the system. At lowest T and n 
particles fluctuate mostly around their mean positions. With in¬ 
creasing T or n the particle diffusion grows. 


We can conclude, that the applicability of the KTHNY 
physics in the phase diagram, cf. Fig. 0}:, is limited to 
the region of thermal melting. In the quantum domain, by 
increasing the system size, we expect to find a first order 
phase transition. Its identification lies in a sharp increase 
of the number of defects, as found for circular and irregu¬ 
lar confinement (admittedly, a similar sharpness is missing 
for the irregular geometry). In contrast, if a transition is 


of the KTHNY-type, the defects are increasing continuously. 
In our case, this crucial distinction is masked by the solid- 
liquid crossover due to finite size effects. In the quantum 
case (which in the thermodynamic limit would correspond 
to a first order transition) we observe a more rapid increase 
of the defect fraction compared to the thermal melting, cf. 

Fig. ini 

To further elaborate possible differences in the melting 
mechanism, we present the snapshots of particle trajecto¬ 
ries in the irregular and elliptical traps during the quantum 
and classical crossover, cf. Fig.jT] In the quantum case, only 
the BC-particle coordinates o are studied to make compar¬ 
ison on equal grounds. It is evident from Fig. [7] that the 
quantum melting sets in primarily by the harmonic fluctu¬ 
ations around equilibrium positions, and the mean squared 
displacements are increased with n. In contrast, the thermal 
melting commences by connecting neighboring particles over 
long distances with coherent displacements making a long 
and tortuous paths of collective motion in certain regions of 
the system. Such a collective motion across the crossover 
also occurs for the circular and elliptic traps, however, the 
paths are less rambling and follow the underlying confine¬ 
ment symmetry. 


G. Distribution of C with classical and quantum 
fluctuations 

Our results from the previous section indicated that the 
quantum crossover from the ’solid’ to ’liquid’ occurs via 
harmonic fluctuations, whereas the corresponding classical 
‘melting’ begins by aligning displacements of certain par¬ 
ticles over long distances while the other particles remain 
confined around their equilibrium positions. We expect to 
illustrate this key difference using the distribution of Linde- 
mann parameter obtained from the statistics over all parti¬ 
cles. This is because a harmonic melting produces a Gaus¬ 
sian distribution of (positive) C centered at zero, but our 
picture of thermal melting would produce its broad distri¬ 
bution with additional peaks (corresponding to the drifting 
particles along the tortuous path of delocalization) and with 
no obvious symmetry. With this motivation, we calculate 
Ci separately for each particle i (f = 1 to fV) by averaging 
it over all the MC configurations. We then represent their 
probability density of displacements as P{C) by construct¬ 
ing (normalized) histograms with these Ci. For our calcu¬ 
lations, Ci is defined as Ci = \/(|ui]y for ith particle, and 

Ui = Ti — . Here is the position of the ith particle in 

the initial configuration (after equilibration). For a justified 
comparison of the results from quantum crossover with the 
corresponding classical ones, the quantum version of Ci is 
evaluated using the BC-coordinates, see Eq. 4. 

Our results for the evolution of P{C^'^) for all the three 
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FIG. 8: Evolution P{C^^) vs n at T = 0.005 is presented where 
BC is used to calculate the distribution to compare it with the 
classical case presented in Fig. [9] In circular case, the distribution 
is narrowly peaked at some high values of n, while in the irregular 
case, P{£^^) is rather broad starting from very low values of 


traps are presented in Fig. [51 In the absence of a good res¬ 
olution of our data, we focus on the the important qual¬ 
itative features of P{C^^) leaving out their details (note 
that a total of IV = 57 particles is not enough for draw¬ 
ing quantitative conclusions). Considering first the irregu¬ 
lar trap (left panels), we notice that the distribution has a 
sharp peak at « 0 for small n (cf. n = 0.06) rep¬ 

resenting the quantum ‘solid’. The peak broadens yielding 
a longer tail with increasing n, however, still a significant 
fraction of particles contribute to fn 0 even beyond the 
crossover, cf. n = 0.125 > nx (note, that for the irregular 
trap nx ~ 0.10(1), cf. Fig. 2d). While most particles tend to 
diffuse in the quantum ‘liquid’ phase, the extent of diffusion 
is not the same for all particles. Particularly the ones near 
the corners and long limbs of the irregular trap, keep ‘rat¬ 
tling’ only around their mean positions contributing weight 
near k. 0. The structure and evolution of P{C^^) for 
irregular trap is qualitatively consistent with our picture of 


FIG. 9: The evolution of classical P{C) with increasing thermal 
fluctuations in the three traps. At large T — 0.045, beyond ther¬ 
mal crossover the distribution is narrow and symmetric in the 
circular case, while in the irregular and elliptical case, it has a 
long-tail like structure. The influence of underlying geometry of 
the trapping potential can be seen from this criteria. 


“melting via harmonic fluctuations”. 

The elliptical and circular traps, on the other hand, pos¬ 
sess spatial symmetry that results in spatial shells of popula¬ 
tion of particles along the symmetry directional^. Thus, for 
the circular and elliptic traps P{C^^) provides also a mea¬ 
sure of the angular shell rotation. Detailed analysis of our 
data demonstrates that the shift of the centroid of P{£^^) 
with n from r; 0 to a finite value, leaving no weight 

to P{C^^ « 0), for circular and elliptic traps is due to the 
loss of the angular order of shells. On the other hand, the 
shape of P{C^^) and its width is contributed by the inter¬ 
shell hops of the particles. The absence of spatial shells for 
the irregular confinement makes the evolution of P{C^^) 
qualitatively different from those in the traps with spatial 
symmetries. 

The evolution of P{C) in Fig. 9 with T in case of thermal 
melting is similar to Fig. 8 insofar as the broad features are 
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concerned (particularly for the elliptic and circular traps), 
though important differences persist. For example, P{C,) is 
already broad for the irregular trap at low T = 0.01. An 
increase of T produces larger fraction of particles whose C 
are distributed over a wide range. Such evolution is in sharp 
contrast with the “harmonic melting” of irregular quantum 
‘solid’ discussed above. Nevertheless, we close this section 
noting that it is the evolution of P{C) that points towards 
differences in the ‘melting’ coming from the geometry of the 
confinement, which is not discerned from most other observ¬ 
ables. 


IV. CONCLUSION 

In conclusion, we have studied the structural order- 
disorder transition for classical and quantum IWM of 
trapped Coulomb-interacting Boltzmann particles. Confine¬ 
ments with or without spatial symmetries leave signatures 
somewhat differently on the crossover from Wigner molecule 
to liquid-like phase. However, the quantum zero point mo¬ 
tion seems to screen the boundary effects in all traps, which 
results in a broadly similar evolution of observables across 
the crossover defining different melting criteria. 

We have introduced several independent melting criteria 
which for different symmetries of the trap point towards a 
unique density interval of the crossover for the quantum melt¬ 
ing, i.e. [nx — An, nx + An] with nx ~ 0.11 being the 
mid-point and An 0.01. This interval via Eq. can be 
translated into the range of bulk densities: 46 < < 69. 

The critical value, Vg « 66.5, found for the quantum melting 
of a 2D Coulomb bulk system of Boltzmannons^ is close to 
the upper bound of the crossover region. A lower crossover 
scale for a trap than the bulk critical value is expected^, be¬ 
cause the loss of translational symmetry in traps reduces the 
ability for the particles to delocalize and the resulting ‘solid’ 
survives up to a smaller r^. 

Our detailed analyses of melting in the irregular confine¬ 
ment allowed us to reconstruct a precise phase diagram which 
shows several phases: a) thermal quasi-classical liquid, b) 
Wigner solid and c) quantum liquid. Here we note that 
the properties of the liquid phase strongly depend on the 
quantum statistics. This will modify the absolute values of 
(nx. An) characterizing the crossover. However, our discus¬ 
sion of the melting criteria introduced for Boltzmann parti¬ 
cles will remain valid also for fermionic and bosonic IWM. 

Of course, boltzmannons are just a model case to capture 
quantum diffraction effects. This neglects quantum statis¬ 
tics effects of fermions or bosons. However, in the crystal 
phase where particles are strongly localized, these effects are 


expected to be smalli. In particular, the suppression of the 
fermi liquid behavior by the strong interactions is well-known 
in the ^He system^. The statistics effects become more pro¬ 
nounced near the solid-liquid crossover line. For ^He systems 
long exchanges of indistinguishable bosons play an important 
role in the stabilizing the superfluid phase^S. In contrast, 
Fermi statistics is known to stabilize the solid phase which 
is directly reflected in the different critical rg-values of crys¬ 
tallization in 2D bulk: Vg ~ 37 for fermionsi^ and Vg Ri 60 
for bosons^. As a consequence, the statistics has an effect 
on the precise location of the crossover line, and should def¬ 
initely influence a location of the crossover region in finite 
systems^. At the same time, we expect that the role played 
by quantum statistics will be reduced due to finite size and 
geometry effects. 

It is indeed reassuring that the phase boundary produced 
by independent criteria match reasonably well within the sta¬ 
tistical errors. We have uncovered a new feature in the phase 
diagram for the irregular trap in the form of the re-entrant 
solidification or the weakening of liquidity upon increasing 
temperature in the quantum regime. Identification of such 
a phenomenon would be an intriguing experimental proposi¬ 
tion. 

The quantum crossover is associated with a sharply in¬ 
creasing defect fraction (disclination) in the region around 
nx- The sharpness in the proliferation of defects is weaker 
in the irregular trap, however, the presence of spatial symme¬ 
tries in case of circular or elliptic confinements strengthens 
such sharpness hinting that the crossover turns into a first 
order transition in the thermodynamic limit. Nevertheless, 
our extrapolation based on a small system can hardly differ¬ 
entiate between a weakly first order and nearly second order 
transitions. They pose serious numerical challenges even for 
a bulk system. In contrast, the corresponding classical tran¬ 
sition is governed by the KTHNY physics. This suggests 
that the quantum melting mechanism should be quite simi¬ 
lar to the one in a harmonic solid in 3D. We hope that our 
work will contribute to pivotal understanding of the melting 
in disordered 2D materials and will stimulate a new body of 
experimental research. 
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